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Abstract 

We study the numerical recovery of Manning's roughness coefficient for the dif- 
fusive wave approximation of the shallow water equation. We describe a conjugate 
gradient method for the numerical inversion. Numerical results for one-dimensional 
model are presented to illustrate the feasibility of the approach. Also we provide 
a proof of the differentiability of the weak form with respect to the coefficient as 
well as the continuity and boundedness of the linearized operator under reasonable 
assumptions using the maximal parabolic regularity theory. 
Keywords: diffusive shallow water equation, parameter identification 

1 Introduction 

The diffusive wave approximation (DSW) of the shallow water equations (SWE) is often 
used to model overland flows such as floods, dam breaks, and flows through vegetated 
areas [20, 12, 8]. The SWE result from the full Navier-Stokes system with the assumption 
that the vertical momentum scales are small relative to those of the horizontal momen- 
tum. This assumption reduces the vertical momentum equation to a hydrostatic pressure 
relation, which is integrated in the vertical direction to arrive at a two-dimensional sys- 
tem known as the SWE. The DSW further simplifies the SWE by assuming that the 
horizontal momentum can be linked to the water height by an empirical formula, such as 
Manning's formula (also known as Gauckler- Manning formula [9]) [4, 21]. The DSW is a 
scalar parabolic equation which resembles nonlinear diffusion. 
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The DSW gives rise to the following initial/boundary value problem for the water 
height u 

du 

— - V • (k (u, Vu) Vu) = / in Q x (0, T] 

u — uq on x {t = 0} (i) 

(fc (u, Vu) Vu) • n = h onT N x(0,T] 

u = g on x (0,T] 

where Q is an open bounded domain in M d (d = 1,2), and Tn and Td are disjoint subsets 
of the boundary T = dQ such that T = Tjv U Tp. The forcing function (e.g., rainfall 
acting as a source or infiltration acting as a sink) / : fi x (0, T] — > R, the initial condition 
Uq : Q — > R, and the Neumann and Dirichlet boundary conditions h : X (0, T] — > R 
and g : x (0, T] — » R are given. The diffusion coefficient Vw) is given by 

c/ |Vtt| |Vm| ' 

where z : :— >■ IR + is a nonnegative time-independent function that represents the bathy- 
metric or topographic measurements available for the region under analysis. The param- 
eters 7 and a satisfy < 7 < 1 and 1 < a < 2. Following Manning's formula [ ], we set 
these parameters to 7 = \ and a — |. The function Cf (or equivalently df — ^) repre- 
sents Manning's roughness coefficient, also known as the friction coefficient. The typical 
values are available in the literature [3, 19]. We refer to [2, 17] for recent mathematical 
analysis and to [17, 5] for efficient numerical algorithms. 

In practice, the Manning coefficient c/ is an empirically derived coefficient, and histor- 
ically it was expected to be constant and a function of the roughness only. It is now widely 
accepted that the values of the coefficients c/ are only constant within some range of flow 
rates, and depends strongly on many factors, including surface roughness, sinuosity and 
flow reach. The presence of multiple influencing factors renders a direct measurement of 
the coefficient values less reliable and the use of a single-valued coefficient also greatly 
constrains the practical utility of the DSW model to faithfully capture important physical 
features of real open channel flows, for which a spatially-varying coefficient is necessary 
due to distinct physical characteristics of different regions. 

In this study, we propose to estimate the distributed Manning coefficient directly from 
water height measurements using inversion techniques, that is, formulating an inverse 
problem for identifying the friction coefficient Cf from measurements of the water-height 
acquired by sensors and infrared imaging. In comparison with direct measurement, the 
proposed approach does not require a knowledge of the physical properties of the overland 
environment, which might be difficult to directly incorporate, and moreover, can naturally 
handle spatially varying coefficients. Therefore, a reliable and efficient estimate of this 
coefficient is expected to greatly broaden the scope of the DSW model and to facilitate 
real-time simulation of the flow, which is of immense significance in a number of applica- 
tions, for example flood prediction and flood hazard assessment. The goal of the present 
study is to propose an inversion algorithm and demonstrate its feasibility on simulation 
data for one-dimensional models. 
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We briefly comment on relevant studies on the inverse problem. Due to its conceived 
practical significance, it has received some attention in the literature [6, 7]. For exam- 
ple, Ding et al [6] estimated the Manning's coefficient in the SWE within the variational 
framework using the limited memory quasi-Newton method, and compared its perfor- 
mance with several other optimization algorithms. However, these works have considered 
only the situation of recovering a few parameters (with a maximum three), instead of esti- 
mating a distributed Manning's coefficient like here. If the number of unknowns is small, 
the ill-posed nature of the problem does not evidence directly. Therefore, the present 
work represents a nontrivial step towards the important task of estimating distributed 
Manning's roughness coefficients. 



2 Linearization of the forward map 

In this section we describe the linearization of the forward map F : df — > u(df), where 
u(df) denotes the solution to system (1). The linearization is required for solving the 
forward problem (with a predictor-corrector method) and the inverse problem (adjoint and 
sensitivity problems, see Section 3). Therefore, its derivation is of independent interest. 
In order to make the presentation accessible, we choose to derive the derivative operator 
informally. A rigorous derivation can be found in Appendix A. 
The bilinear form of problem (1) is 



B(u, w) 



u t wdx + / k(u, Vm)Vm ■ Vwdx 
n Jn 

u t , w) + (k(u, Vu) Vu, Vw) , 



and the linear form is 



£{w) 



fwdx + / hwds = (f, w) + (h, w)r r 



The weak formulation of the problem reads: For almost all t G (0,T], find u with the 
given Dirichlet boundary condition and initial data u(0) = Uq such that 

B(u,w)=£(w) VweV, 



where V is an appropriate function space [17]. 

We shall seek the Gateaux derivative of the bilinear form B at u, that is, ^B(u + 
ev,w)\ e= o- We aim at deriving an explicit formula to facilitate further developments. We 
proceed as follows. It follows from the product rule for differentiation that 



dB(u + ev, w) 



e=0 



d_ 
de 



, , \(u+ ev) — z] a „ . „ 

(u t + ev t , w) + d f \h {Vu + eVv) , Vw 

\Vu + evvl ' 



e=0 



( u _ z ) a 

(v u w) + ( d f \_ ,iLy v, Vw) + I + II, 



Vu 



1-7 
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where the terms / and II are respectively given by 



, d\u + ev - z] a Vu + eVv „ 
d f — * — 7^ ^rr, — ,Vw 



|Vw + eVw| 1 -i" 



df a 



u 



>a-l 



e=0 



IVu 



1-7 



v Vu, Vw ) , 



and 



II — ( df[u + ev — z 



^{Vu + eVv^' 1 



de 



(Vu + eVv) , Vw 



e=0 



(d f (u-z) a ( 7 -l)|V« 



7 2 • Vv Vu, Vw 
Vw 



df (7 - 1) 337 Vu ■ Vv Vu, Vw 

\Vu\ 7 



Here the second line follows from the relation |Vu| = VVu • Vu = (Vu-Vu)^ that implies 



d\Vu + eVv\ 



de 



e=0 



\ ((Vu + eVv) ■ (Vu + eVv)) * 2(Vu + eVv) ■ Vv\ e=0 = ^ V . 



Consequently, by combining all these identities, we arrive at the following formula 



dB (u + ev , w) 



de 



e=0 



= (v t ,w)+ (df y ' Vv,Vw 



+ ( dfd 



[u — z 



,a-l 



|Vd 1 -T 



v Vu, Vw 



1 , / — z) a Vu „ Vu T „ 

+ d/(7-l) ;_ ,/ -Vv t^,Vw 



Wul 1 -^ IV 



u\ 



Vu 



(v t , w) + (k(u, Vu) (I — (1 — 7) fj <g) fj) • Vv, Vw) 



+ k(u,Vu) 



a 



(u-z) 



v Vu, Vw 



where / is the identity operator and the vector field fj = is the normalized gradient 
vector field. The matrix-valued function fj ® fj represents a projection operator onto the 
gradient direction fj. Hence, the structure of the second term indicates that, for the 
linearized problem, the diffusion along the gradient direction is attenuated by 1 — 7, 
whereas the tangential component is not affected. To simplify notation we denote this 
attenuated diffusion tensor as 



k vv (u, Vu) = k(u, Vu) (I — (1 — 7) fj <g) fj) . 

Meanwhile, the linearized problem has a convection term (the third term), as a conse- 
quence of the nonlinear term involving u. These structural terms relate to the underlying 
physics of the model. 



4 



It follows directly from the definition of the Gateaux derivative, i.e., which is de- 
noted by v = u'(df)d G V and characterizes the perturbation of u(df) caused by a small 
perturbation of the coefficient df in the direction d that it (in weak formulation) satisfies 

(v t , w) + (k m (u, V«) ■ Vv, Vw + / ' y w V«, Vwi = - d- nbrVu, Vw 

and the initial condition is v (0) = 0, since the initial data is not affected by a perturbation 
of the friction coefficient. 



3 Inversion algorithm 

Now we turn to the inverse problem of reconstructing the coefficient df from the measure- 
ments of water heights. As a general rule, the inverse problem is ill-posed in the sense 
that small perturbations in the data can lead to large changes in the solution. Hence we 
adopt a regularization strategy by incorporating a penalty term into the cost functional, 
following the pioneering idea of Tikhonov [ ] . More precisely, we consider the following 
penalized misfit functional 

J(d f ) = - f f (u(d f ) - g) 2 dxdt + - f \Vd f \ 2 dx, 
2 Jo Jn 2 

where the scalar 5 is the regularization parameter, and g denotes the noisy measurements 
of the water height u(df). With minor modifications, the algorithm discussed below can 
also be applied to other measurements, for example, water height on the boundary or 
scattered in the domain. The term ||Vd/||| 2 (-m enforces smoothness on the sought-for 
coefficient, and thereby restores the numerical stability necessary for practical computa- 
tions. To numerically minimize the functional, we adopt the conjugate gradient method. 
The method is of gradient descent type, and it only requires evaluating the gradient of 
the functional J(df) at each step. We note that the conjugate gradient method has been 
successfully applied to a wide variety of practical inverse problems, such as in heat transfer 
and mechanics; see for example, [1, 16] and references therein for details. 

To derive a computationally efficient gradient formula, we first note that, given a 
(descent) direction d, the misfit term in the functional J can be approximated using a 
Taylor expansion and ignoring higher order terms. 



1 ' T 



- / / ( u (df + d) — gfdxdt — - f f (u(df) — g) 2 dxdt = 
2 io Jn % Jo Jn 

u(df + d) — u{df)) (u(df + d) — g + u(df) — g) dxdt 
u'(df)d (u(df) — g) dxdt. 



2 



o Jn 



'o Jn 

The approximation is reasonable if the magnitude of the direction d is small. 
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The last formula can be further simplified with the help of the adjoint problem for p, 
which in weak form reads 

/ \ n i t-, \ t-, t-, \ ( o>k(u,Vu) \ , 

(-pt, w) + {k m (u, Vu) ■ Vp, Vw) + I j~ ^ Vu ■ Vp, w I = {u(df) - g, w) 

together with the terminal condition p(T) = 0. Recall the weak formulation of the 
sensitivity problem v — u'(df)d, that is, 

/ x /, , „ s „ „ s fak(u,Vu) „ „ \ (An— z) a „ „ 
(v t , w) + (k m (u, Vu) ■ Vv, Vic) + f ^ ^ v Vu, Vw\ = - f d K ' Vu, Vw 

together with the initial condition v(0) = 0. Upon setting the test function w = u'(df) d 
and w = p in the weak formulations for p and u'{df) d, respectively, we arrive at 

/ (u(df) — g,u'(df) d) dt = — / / d . Vp-Vudxdt— I —(p,u'(df)d)dt 
Jo Jo Ja 7 Jo dt 



T 



JO, 



d K — . z ' Vp ■ Vudxdt, 
'Vm 1 " 7 



where the last identity follows from the initial condition for u'(df) d and terminal condition 
for p. This relation yields the following concise gradient formula of the functional J{df) 



j'(df) = - / ) r^-Vp ■ Vudt - 6Ad f . 



We note that this gradient J'{df) is inappropriate for updating the coefficient df directly 
due to its lack of desired regularity. The consistent gradient of the functional with respect 
to H\to), denoted by J' s (d f ), can be calculated as 

-AJ' s (d f ) + J' s (d f ) = J'(d f ) 

with a homogeneous Neumann boundary condition. 

Now we can give a complete description of the conjugate gradient method summarized 
in Algorithm 1. In the algorithm, one has the freedom to choose the conjugate coefficient 
(3^ and the step size 9k- There are several viable choices of the conjugate coefficient [11]. 
One popular choice is suggested by Fletcher-Reeves, which reads 

ll^(^/)ll|2(Q) 

Pk-l ~ 



iwnin^) 

with the convention f3 = 0, and then update the conjugate direction dk with 



d k = J' s {d)) + /3 k -id 



k-i- 



Generally, the step size selection is of crucial importance for the performance of the 
algorithm. We have opted for the following simple rule. By means of a Taylor expansion 
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Algorithm 1 Conjugate gradient method. 

1: Set k = and choose initial guess d®. 
2: repeat 

3: Solve direct problem with df = d k , and determine residual r k = u(d k ) — g. 
4: Solve adjoint problem with right hand side r k . 

5: Calculate gradient J' s (d k ), conjugate coefficient (3 k , and direction d k . 

6: Solve the sensitivity problem with direction d = d k . 

7: Compute step length 9 k in conjugate direction d k . 

8: Update coefficient d k = d k — 9 k d k . 

9: Increase k by one. 
10: until A stopping criterion is satisfied. 
11: Output approximation df 



of the objective function J(d k — 9d k ), with the forward solution u(d k — 9d k ) linearized 
around d k , we arrive at the following approximate formula for determining an appropriate 
step size k 

e _ (r fc , u'(d k f )d k ) L 2 {0tT . L 2 {n)) + 5(Vd k f , Vd k ) L 2 {n) 

\W(df)dk\\L2(o,T;L2(ci)) + <5||V4||i 2(n) 

where r k = u(d k ) — g denotes the misfit (residual). The step size 9 k is determined to 
enforce a reduction the functional value, that is, J [d k — 9 k J' s (d k )) < J(d k ). Our ex- 
perience with other inverse problems indicates that this choice works reasonably well 
in practice [16]. Advanced step size selection rules, such as, Barzilai-Borwein rule with 
backtracking, maybe also be adopted to further enhance the performance. The algorithm 
terminates if the selected step size falls below 1.0 x 10~ 3 . Overall, each step of the iter- 
ation invokes three forward solves: the (nonlinear) forward solve for computing the map 
u(df), the (linear) adjoint solve for calculating the adjoint p(df) and consequently the 
gradient J'{df) and the (linear) sensitivity solve for selecting the step size 9. The extra 
computational effort for computing the smoothed gradient J' s {df) is marginal compared 
with other steps due to its simple structure. 

4 Numerical experiments and discussions 

Here we present some numerical results for one-dimensional examples to illustrate the 
feasibility of the proposed inversion technique. The forward problem is discretized using 
piecewise linear finite elements in space and the generalized-a method in time (detailed in 
Appendix B). The adjoint and sensitivity problems are both solved with the generalized-a 
method. 

The spatial domain £7 = [—2, 2], and the mesh size h is |. The time interval is [0, |] , 
and the time step size is ^. This mesh was used for both generating the exact data 
and used in the inversion step (i.e., adjoint problem and sensitivity problem). We note 
that we also experimented with using finer mesh for generating the exact data, and the 
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Table 1: Numerical results (error e) for different noise levels. 



0% 0.5% 1% 2% 



Example 1 5.94e-3 2.00e-2 2.90e-2 4.52e-2 
Example 2 4.49e-2 6.41e-2 7.49e-2 9.99e-2 
Example 3 4.05e-2 5.15e-2 5.94e-2 9.42e-2 



reconstructions are identical. Also both the forward solution u(df) and the coefficient df 
are represented in this mesh. The initial guess for the coefficient is df = 1. The noisy 
data g are generated pointwise as 



where e is the relative noise level, and the random variable (noise) ( follows a standard 
Gaussian distribution. The choice of the regularization parameter 5 is crucial in any 
regularization strategies [18]. There have been intensive studies on its appropriate choice 
which have led to systematical and rigorous rules for choosing an appropriate value, 
see [15, 13] for recent progress. However, in this preliminary study, we have opted for the 
conventional trial-and-error approach. 

We consider three examples: one with a smooth coefficient, and two with a discontin- 
uous coefficient. First, we consider the recovery of a continuous coefficient. 

Example 1. The forward problem has a homogeneous Neumann boundary condition, and 
the initial condition w is Uq — — \x + |. The exact coefficient is d\(x) = 1 + jq(x 2 — 4) 2 . 

Figure 1(a) and Table 1 show the numerical results for Example 1, where e is the 
relative error of an approximation df, defined as e = \\df — dj||L 2 (n)/IM/||L 2 (n)- The 
reconstructions are in reasonable agreement with the exact coefficient df for up to 2% 
noise in the data. Hence the proposed method is stable and accurate. We note that 
the approximation near the boundary seems less accurate compared to other regions. 
The error e decreases as the noise level e decreases to zero, see also Table 1. Overall, 
the convergence of the inversion algorithm is rather steady, see Figures 1(b) and (c). 
While the functional value J(d k f) decreases monotonically as the iteration proceeds, the 
convergence of the error e exhibits a clear valley, indicating that a premature termination 
of the algorithm might result in sub-optimal reconstructions. 

Then we consider the recovery of a discontinuous coefficient. 

Example 2. The boundary condition and initial condition of the problem are identical 



characteristic function. 

Figure 2(a) and Table 1 present the numerical results for Example 2. The convergence 
of the result with respect to the noise level e is again clearly observed. The reconstructions 
capture the overall shape of the true solution. However, the discontinuities are not well 
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(a) reconstructions 



20 40 60 80 100 
k 

(b) functional value J{dh) 



20 40 60 80 100 

k 



c error e 



Figure 1: Numerical results for Example 1. Here the convergence of Algorithm 1 is for 
e = 2% noise. 
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(a) reconstructions 
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(b) functional value J{dt) 




c error e 



Figure 2: Numerical results for Example 2. Here the convergence of Algorithm 1 is for 
e = 2% noise. 



resolved, even for exact data, and consequently the results are less accurate compared 
with those for Example 1. This is attributed to the presence of discontinuities in the 
sought-for solution dj, which cannot be accurately approximated using the smoothness 
penalty |Vd/|^ 2 (nv Discontinuity preserving penalties, such as, total variation, might 
be employed to improve the resolution. Nonetheless, the conjugate gradient algorithm 
remains fairly steady, see Figures 2(b) and (c). 

A last example considers the recovery of a more complex coefficient profile. 

Example 3. The boundary condition and initial condition of the problem are identical 
with those in Example 1. The exact coefficient d^(x) is given by dj = 1— \x\-Z -^]+§X[^ aj- 

Here the true solution has more refined details, and hence the spatial mesh size h is 
accordingly refined to | for a better resolution. The results for Example 3 are shown 
in Figure 3 and Table 1. The convergence of the numerical reconstruction with respect 
to the noise level is again observed, see Table 1. The observations for the previous 
example remain valid: the numerical reconstructions roughly capture the profile of the 
true solution, but fail to resolve accurately the discontinuities, and the algorithm converges 
steadily and reasonably quick. 
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x 10 20 30 ' 10 20 30 

k k 

(a) reconstructions (b) functional value J(dJf) (c) error e 

Figure 3: Numerical results for Example 3. Here the convergence of Algorithm 1 is for 
e = 2% noise. 

5 Concluding remarks 

We have presented an inversion technique for estimating the Manning's coefficient in the 
diffusive wave approximation of the shallow water equations. The results show that the 
proposed approach is capable of yielding an accurate and stable estimate in the presence 
of noise. We have also detailed a careful study of the properties of the forward map, 
in particular, we discuss its continuity and differentiability based on maximal regularity 
theory for parabolic problems. The mathematical analysis, such as, convergence and 
convergence rates, of such an inversion technique remains to be investigated. Also the 
evaluation of the method on real data is of significant interest. 
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A Properties of the forward map 

In this part, we briefly discuss the continuity and differentiability of the forward map 
F : L°° — > L 2 (0, T; df i— > u(df) based on maximal regularity theory for parabolic 

problems [10]. The conditions in Theorem A.l impose a certain regularity constraint 
on the coefficient df as well as on the boundary and initial conditions. Such mapping 
properties are essential for analyzing commonly used regularization schemes, for example, 
Tikhonov regularization and Landweber iteration for solving the inverse problem, and for 
establishing the convergence of numerical algorithms. 
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We first show the Lipschitz continuity of the forward map F. 

Theorem A. 1. Assume that \\u(df)\\L°° and\\Wu(df)\\Loo are uniformly bounded, and that 
df, \Vu(df)\ and (u(df)-z) are strictly positive, and further the gradient \V(tu(df) + (1 — 
t)u(df))\ is strictly positive for all t G (0, 1) and df, df in the admissible set A. Then if 7 
is sufficiently close to unity, the mapping F : L°° — > L 2 (0, T; given by df >->■ 

zs Lipschitz continuous on A. 

Proof. We denote by fc (m, Vm; df) = df |y~|i- 7 and m = u(df), and let t> = m — m. We 
denote the bilinear form parametrized by df as 

B(u, w; df) = (u t , w) + (k(u, Vm; df) Vm, Vw) , 

By subtracting the bilinear forms B(u, w; df) = £(w) and B(u, w, df) = £(w) and choosing 
w — v, we arrive at 

= (ut — Ut, w) + (k{u, Vm; df)Vu — k(u, Vm; d/) Vm, ViuJ , 

which by virtue of the assumptions on u and Vm can be rearranged into 

2^IMIl 2 + C k\\Vv\\ 2 L 2 < - (k(u, Vm; d f - d f )Vu, Vu) 

- ( (fc(u, Vm; d f ) - k(u, Vm; d f )j Vm, Vu) := J + II, 

where Ck is the coercivity constant for the bilinear form B(-,-). Using Cauchy-Schwarz 
inequality and Young's inequality, the first summand / on the right hand side can be 
estimated as follows 

I<C(e 1 )\\df-df\\ 2 LOO + e 1 \\Vv\\l 2 . 
Meanwhile, we split the nonlinear term in the bracket in the second summand II into 



k(u, Vm; df) — k(u, Vm; df) = df 



(m-^hivmI 


1-7 _ 




Vm 






Vm 


i-^VM 


1- 


-7 





+ 



|Vm| 1 -t 



(2) 

Now the mean value theorem gives 

(m — z) a — (m — z) a = a (M — z) a ~ l v, (3) 
where M is an element between u and u, and also by means of the Taylor expansion 

IVm) 1 " 7 - IVm] 1 " 7 = (1 - 7 )v • Vv, (4) 



and the function 

" x V(u-sv) 



o |V(m- sv)\ 1+ "< 



ds, 
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which by assumption is bounded in L°°. Consequently by Young's inequality, we get 

II < (1 -7) \\k{u, V«, d/)|U- HvlUoo 1| VtiHIoo || VwHia + C(eJ>||£ a + f ||Vu||f 2 ) 
<C ((1 - 7) + e 2 ) UVullla + Ce^lMI^. 

Since 7 is close to unity and for sufficiently small ei, e 2 , fi '■— C ((1 — 7) + e\ + e 2 ) < Ck, 
we obtain 

§&Nl£a + {C K - n) \\Vv\\ 2 L 2 < Ce^lMlla + C||d/ - rf/lli- 
Now an application of Grownwall's inequality leads to 



Jo 



upon noting the condition u(0) = u{0). □ 

Our next result improves the regularity of the map in Theorem A.l by invoking 
Groger's maximal regularity theory [10], which is needed for the differentiability. 

Theorem A. 2. Let the assumptions in Theorem A.l be fulfilled. Then the mapping 
F : L°° — > L 2 {0, T; W 1,P {Q)), df 1— > u{df) is Lipschitz continuous for some p G (2, 00). 

Proof. As before, we denote by 

, , _ , . , {u — z) a 
klu, Vm; df) = df — M — 
v 11 1 IVmI 1 ^ 

and u = u{df), and let v — u — u. Then v solves 

v t + Av = f 

with 

Av = -V • ( (k{u, Vm; If) - k{u, Vm; d f ) \ Vtt + k{u, Vu; d f )Vv^j 

and / = V- (k{u, Vu; d f - d f )Vuj . Clearly, / e L p {0, T; {W 1 *)') for d f -d f e LP because 

the remaining terms are uniformly bounded in L°°. To apply Groger's theorem [10, 
Theorem 2.1], we only need to show the coercivity and boundedness of the operator A 
defined above. By using the Taylor expansions (3) and (4) in the splitting (2), we can 
rearrange the differential A into 

{Av, w) = (k{u, Vu; df) (1 - 7) | Vtt| 7_1 v • VvVu, Vu?j 

- (d f a {u - zf 1 |Vu| 7 ~Wu, V«j) + (k{u, V«; d f )Vv, W) . 

In view of the strict positivity of the term k{u, Vit; df), that the parameter 7 is close to 
one and that the quantities u, Vu etc. are uniformly bounded, we deduce that 

{Av,v) > c A \\Vv\\ L 2 - C A \\v\\ L 2. 



12 



for some constants ca, Ca > 0. Hence, the associated matrix-valued coefficient in the 
differential operator is pointwise bounded from below and above away from zero. The 
continuity of the operator follows similarly. Consequently, an application of Groger's 
theorem [10] directly yields the desired estimate f Q \\v(s)\\^i tP ds < C||d||_£,cx> for some 
j9G(2,oo). □ 

Remark A.l. The exponent p e (2, oo) in Theorem A. 2 depends on the spatial dimension, 
the pointwise upper and lower bounds of the conductivity k(u, Vu; df) and the smoothness 
of the domain Q; see [10] for details. 

Next we show the boundedness of the linearized map. 

Theorem A. 3. Let the assumptions in Theorem A.l be fulfilled, and the linear map 
F' : L°° — > L 2 (0,T; H 1 ^)) be defined by e£ i-» v, with v given by 

(v t , w) + (A; (it, Vu) [I — (1 — 7) fj <g) fj] • Vf , Vw) 

t(u — z) a_1 t> _ \ / ,(u — z) a _ _ 
J IVmI 1 -^ J V IVuI 1 -) 1 

with the initial condition v(0) = 0. Then the linear map F' is bounded. 
Proof. Insert w := v to get 

§3tlM&(n)+ Vu; d f )Vv, [I - (1 - 7) fj ® fj]Vv) 

d f a^——t Vu, Vv - d „ M 7 Vit, Vw := J + i7. 

Using Cauchy-Schwarz inequality and Young's inequality, the term / can be bounded by 
I < ^eOlld/Hiooll (u - z)*- 1 |V^n|ioo|b||| 2 + ei||V^||i 2 , 

where 6\ > is arbitrary. Similarly, the term II can be bounded by: for any t 2 > 
// < CCeaJUdllioolld/lliooll (u - z) a |Vu| 7 ||i 2 + e 2 ||Vv||ifl. 

Recall that 7 is strictly less than unity, and hence I — (1 — 7) rj ® rj is strictly positive 
definite, and the diffusion coefficient k(u, Vu;d/) is strictly positive (independent of df). 
Therefore, these estimates altogether give 

|WI! 2 (n) + l|v^||! 2 <c(|H|i 2 + ||C-)- 

Applying Gronwall's inequality and noting that v(0) = 0, the desired assertion follows. □ 

Remark A. 2. The condition that the parameter 7 is close to 1 is not required in The- 
orem A. 3. A direct application of Groger's theorem indicates that the map F' : L°° 1— > 
L p (Q,T;W 1 ' p (£l)) is also bounded. 
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Finally, we show the Frechet differentiability of the forward map. 

Theorem A. 4. Let the assumptions in Theorem A.l be fulfilled, and the bounded linear 
map u'(df)d be defined in Theorem A. 3. Then u'(df)d is the Frechet derivative of the map 
df — > u{df), i.e., 

lim IKj/ + d) - u{d f ) - u'(d f )d\\ L 2 (0jT . H i m _ 



Proof. We denote by k(u, Vm; df) = df ^ u ^l J and df = df + d, u = u(df), u = u(df) and 
u = u'(df)d and let v = u — u, w = u — u — u. We also denote D{u) — (1 — 7) 77 <2) 77. Then 
it directly follows from the weak formulations for u, u and u that 

(wt, w) + ( k(u, Vw; df) Vw — k(u, Vw; df)Vu, Vw 



(k{u, Vu; d f ) (I - D(u)) Vw, Vw) + [ d f a v iy^jT^ Vm > Vu; 

which upon rearrangement and noting the assumptions on u and Vw yields 
(w t , w) + (k(u, Vw; df)Vw, Vw) = ((k(u, Vw; df) — k(u, Vw; df))Vu, Vw) —(k(u, Vw; <i)Vw, Vw>) 



/ (w — £) a-1 w \ 
{k(u, Vw; d f )D(u)Vu, Vw) + ( d f a K ' Vw, Vw J . 

Ill V v ' 



iv 

It suffices to estimate the four terms on the right hand side. First, by means of Cauchy- 
Schwarz inequality and Young's inequality, the term II can be estimated by 

|//| <C||rf|| L oo||Vn|| L 2||Vw|| L 2, 

To bound the first term /, we further split it into 

T /, , „ ^ A\Vu\ 1 ^ -\Vu\ 1 ^)^^ \ f 1 (u-z) a -(u-z) a ^^\ 

1= (k(u,Vu;d f )^ ' |W|1 L 7 ' L Vu,Vw) + (d/ '-Vu^wj :=V + VI. 

Now we employ the Taylor expansion 

IVmI 1 " 7 = IVwI 1 - 7 + (1 - 7)|Vwp- 1 Vw • Vv + KVv 2 
with the matrix-valued function K given by 

k = - f\i - 1) ((1 - ^mmmr- 3 + a - ^imr^i) * 
jo 
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and 4>{t) = V(u + tv). With the help of this expansion, we derive that 

V-III= (k(u, V«; d f )D(u)Vw, Vw) + (1 - 7 ) (k(u, V«; d f ) ( ^U^, ~ 



+ (I-7) ( fc(^VM) 7 ,V W j + f ^KVn;^)-— V«,V«; 



v v — 

VIII IX 



Next we estimate the terms on the right-hand side one by one. First, let p be the exponent 
from theorem A. 2 and choose q > 2 such that - + - = \. Then by the uniform L°° 
boundness of u and Vii (also u, Vtt etc.) 

VII = (k(u, Vu; d f ) V ^^ 2 V | Vw) 7 ' 1 (\Vu\ 1 -^Vv + VmQVuI 1 - 7 - |Vm|^ 7 )) , Vwj 

< C||Vw|| LP ||Vt;|| L9 ||Vw|| L 2 +C , |||Vw|(|Vn| 1 ^- IVmI 1 " 7 )!!^!^^!!^ 

< C||Vv||lp||Vt;||l«||Vu'|| L 2 + C|||V?;| 2 ||i,2||VHU 2 

< C\\Vv\\ LP \\Vv\\ Lq \\Vw\\ L 2 < C\\Vv\\}% d \\Vw\\ L 2, 

where in the third line we have utilized the expansion (4), and the last line follows from 
the fact that either ||Vt> < C||Vf \\lp holds for q < p or ||Vt> < C||Vf \\ S LP holds for 
some < 5 < 1 due to the L°°-boundedness of Vu and Vu. Similarly, the terms VIII 
and IX can be bounded by 

VIII <C\\d\\ L oo\\Wv\\ L 2\\Ww\\ L 2 and IX <C\\Vvf£\\Vw\\v. 

Next we combine the terms VI and IV. To this end, we employ the Taylor expansion 

(u - z) a = {u- zf + a{u - z) a - l v + \a{a - l){u - z) a ~ 2 v 2 

with u being some function pointwise between u and u we can estimate. With the help 
of this identity, we arrive at the following splitting 

VI + IV = — ( d /a <"-:>;:> <■■ V») + ( d,a(u - zT-H ( ^ - ^) , V W 



X 



( (u-z) a ~ l \ (~ (n - z) a ~ 2 \ 
- da K tiVii, Vm; - | d,a(a - 1) ' n v 2 Vu, Vw 

V v ' S v ' 

XI XII 

Consequently, by the uniform boundedness of the quantities u, Vu (and u, Vu etc.) and 
Sobolev embedding theorem, we have 

X < C\\v\\$ iP \\Vw\\i?, xi < C\\d\\ L oo\\v\\ m , P \\Vw\\ L 2, XII < C\\Vv\\$ tP \\Vw\\v. 
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These estimates, Young's inequality and that 7 is close to unity (hence D{u) can be made 
arbitrarily small for 7 close to unity) yield 

+ iHVHIis < c (\\d\\U\v\\ 2 whP + IMISX + IIHIIO ■ 

Finally, an application of Gronwall's inequality and Theorem A. 2 lead to 

\M\l*+ [ \\Vw\\ 2 L 2ds < c\\d\\ 2 L + J s 



upon noting the initial condition w(0) = 0. This concludes the proof. □ 

Remark A. 3. An inspection of the proof indicates that the assumptions on the solution 
u(df) and gradient can be greatly relaxed if the parameter 7 = 1. The latter case is 
analogous to the porous media equation, and thus the results are of independent interest. 



B Generalized-ce method 

In this appendix, we describe the generalized-a method. Note that for the full discretiza- 
tion of the forward problem, each time step involves solving a highly nonlinear (and 
possibly also stiff) system. Hence a careful treatment of the time stepping is required. 
To this end, we employ the so-called generalized-a method together with a predictor- 
corrector method [14, 5]. For a first-order system, the method can be stated as follows: 
given (u n ,u n ), find (u n+1 ,u n+ i,u n+af ,u n+am ) such that 

u n+a; = u n + a f (u n+1 - u n ), 

1 u n+a m = U n + OC m \U n -\-X — U n ) , 
u n+ i = u n + At((l - 7)m„ + 7M n+ i), 

where At = t n+ i — t n is the time step size, a/, a m and 7 are real valued parameters of the 
method, and R(u n+af , u n+am ) denotes the (discrete) residual of the nonlinear system. For 
a linear model problem, unconditional stability of the scheme is attained if a m > af > |, 
and a second-order accuracy can be achieved with the choice 7 = | + a m — af [11]. The 
method can be succinctly parameterized by the spectral radius into a one-parameter 
family. Then the parameters a m , af and 7 can be expressed as [14] 

1 3 - Pop _ 1 

a/ "i + Poo' am ~2(i + Poo y 7 "i + Poo - 

A complete description of the generalized-a method is given in Algorithm 2. It is of 
predictor/corrector type with correctors computed by a Newton method, where the su- 
perscript indices indicate the corrector steps within the loop. In our implementation, 
we have set = 0.1, and the tolerance e in the stopping criterion to 1.0 x 10 -6 and 
the maximum number of iterations (Maxlter) to 20. The major computational effort of 
Algorithm 2 lies in calculating the Jacobian matrix K n % ^ +1 for the Newton system, i.e., step 
5. For large-scale problems, iterative solvers, e.g., GMRES or BiCGstab, which requires 
only matrix- vector multiplication, are preferable [14]. 
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Algorithm 2 Gcneralized-a method. 



1: Compute predictor u^li = u n and u^li = -h^u n , and set i — 0. 
2: Set the initial guess of u^l af and i4i+a m as 

«i°| a/ =u n + a f (u% - M n ) and u^l am =u n + « m («i+i ~ «n)- 



3: while % < Maxlter do 

4: Evaluate the Newton residual R^+i = R( u n+a f ^n+a m )- 
5: Calculate the Jacobian 

Kji = T^T + am[a/7At] 



6: Solve Newton system for the corrector Att^+i from iY^^Aw^^ = — -R^+i- 
7: Update the solutions 1^+^ and by 

(i+i) _ (i) A (i) 

u n+a f — u n+l ^ ^"ra+l? 

«i+«L = C 1 - 7 _1 am)«S-i + a m [7Aia / ] _1 (u^ ) / - u n ). 

8: Check the stopping criterion: if H-R^H < stop iteration. 

9: Increase index % — % + 1. 

10: end while 

11: Output the solutions u n+ \ and ii n+1 by 

_ . — 1/ (Maxlter) \ i _ . _ 1 /. (Maxlter) . \ 

u n +\ — u n + a j yu n+af — u n ) ana u n+ i — u n + a m [u n+ctm — u n ). 
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